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1.  INTRODUCTION 


In  order  to  produce  stable  solutions  using  the  MACH  code  \Tautz,  et.  al.,  1987]  for  the  Multiple 
Interactions  in  Plasma  (MIP)  chamber  experiment  simulations  [Morgan,  et.  al.,  1989;  Chan,  et.  al., 
1989],  it  was  found  necessary  to  include  a  trapped  electron  model.  The  physical  assumptions  and  the 
numerical  techniques  used  to  implement  this  model  are  presented  in  this  report. 

We  first  describe  the  nature  of  the  instability  that  is  encountered  in  the  chamber  simulations  if  the 
trapped  electrons  are  absent.  The  source  of  this  instability  is  illustrated  schematically  in  Figure  1 .  The 
ion  beam  enters  the  chamber  at  the  throat  and  creates  a  region  of  positive  charge.  This  positive  space 
charge  attracts  electrons  from  the  neutralizing  source  near  the  throat  entrance,  causing  the  electron 
density  to  decrease  (for  strictly  one-dimensional  motion,  the  density  falls  off  inversely  as  the  square  root 
of  the  potential).  On  the  other  hand,  the  ions  are  repelled  by  the  positive  charge  and  tend  to  build  up 
at  the  turn-around  point,  thereby  increasing  the  ion  density.  These  two  effects  produce  a  region  of 
increasing  net  positive  charge  which  does  not  stabilize  until  a  wall  of  positive  potential  builds  up  to  the 
stagnation  energy  of  the  ions.  At  this  point,  the  potential  would  shut  off  the  beam  ion  current.  Since 
this  phenomena  is  not  observed  in  the  chamber,  it  must  be  an  unphysical  effect  due  to  the  neglect  of 
some  essential  physics  component. 

The  MACH  program  normally  proceeds  to  find  solutions  by  Poisson-Vlasov  iteration,  i.e.  the  plasma 
is  assumed  to  be  collisionless.  ('Hie  code  was  designed  for  space  simulations,  where  this  can  be  a  good 
approximation.)  However,  when  applied  to  the  MIP  chamber  simulations,  the  above  instability  problem 
was  encountered.  To  achieve  reasonable,  stable  solutions,  we  have  included  collisional  effects  for  the 
electrons.  That  is,  we  assume  some  electrons  will  scatter  into  the  potential  energy  well  that  is  created 
by  the  streaming  ions,  and  quench  the  unstable  region  of  positive  charge.  We  further  assume  that  the 
electrons  establish  an  isotropic  equilibrium  within  the  well,  with  temperature  equal  to  that  of  the  Vlasov 
(untrapped)  electrons.  The  details  of  this  trapped  electron  model  are  described  below.  In  section  2,  we 
outline  the  "usual"  method  of  obtaining  density  in  MACH  by  means  of  inside-out  trajectories  and  velocity 
space  quadratures.  This  serves  to  introduce  the  basic  equations  and  notation.  In  section  3,  we  describe 
the  modifications  which  enable  us  to  calculate  the  density  of  the  trapped  electrons.  An  approximate 
expression  for  the  untrapped  population  is  also  given.  Density  normalization  of  the  untrapped 
component  is  determined  by  equating  ion  and  electron  fluxes  in  the  source  region.  'Hie  normalization 
constant  for  the  trapped  component  is  essentially  a  free  parameter  of  the  model.  In  section  4,  we  give 
some  numerical  results  for  the  model  for  MIP  experiment  simulation.  Hie  conclusions  are  presented 
in  section  5. 
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Figure  1.  Schematic  of  instability  in  electrostatic  potential 
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2.  THE  INSIDE-OUT  METHOD 


In  the  MACH  program,  the  dimensionless  particle  number  density  n  =  N/N^  is  calculated  by  a  velocity 
space  integral  over  the  distribution  function  /(tT)  at  every  point  r  of  the  spatial  grid 


”  =  JIJ 

=  JJJ  f{v,Q,^)v^dvdQ 


where  dQ  \s  the  solid  angle  element 

dO  =  sin(6)ded4> 

We  now  make  a  transformation  from  velocity  v  to  total  energy  E,  given  by 

E  =  imt>2  +  qV 
2 


where  V  =  V{i)  is  the  local  electrostatic  potential.  By  differentiating  this  equation,  we  get 

dE  =  mvdv 


so  that 


v^dv  =  vfmdv  =  2''^m-^{E-qVy'^ 
and  the  density  integral  goes  into 

n  =  2^^m--'^jjjf{E,e,^){E-qV)^^dEdQ 

By  Liouville’s  theorem,  the  distribution  function  is  constant  along  a  particle  trajectory,  so  that  by 
following  time-reversed  tracks  to  infinity,  we  can  obtain  f(E,  8,  <p)  =  /„  If  we  now  specialize  to  a 
drifting  Maxwellian  distribution  at  infinity,  we  can  write  /„  as 

L  =  (m/(2nikT))^exp(-Jlf) 


where  the  exponent  has  the  form 
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X  =  -  m{v^-v.flkT 
2 


=  e  +  A/^-2e''^-M-cos(0,) 


and 

=  track  velocity  at  infinity 
6.  =  track  polar  angle  at  infinity 

source  flow  velocity  at  infinity 

e  =  EIkT  =  i  m(vJ^/kT 
2 

M  -  t;,/(2JfcT/m)’^  =  Mach  number 
Using  these  variables,  we  can  write  n  as 

n  =  —  jr"^JJJexp(-X)(e- 


where  x  =  is  the  dimensionless  potential  energy.  This  is  the  expression  used  in  MACH  to 

calculate  density  by  means  of  quadratures.  In  order  to  evaluate  the  integrand  of  the  triple  integral, 
particles  are  tracked  backwards  in  time  to  the  source  to  obtain  0^  and  hence  /„  (any  particles  hitting 
objects  contribute  zero  density).  This  tracking  must  be  done  for  each  value  of  e,  6,  (p  with  sufficient 
precision  to  resolve  the  distribution  function  over  all  energies  and  4^  solid  angle.  (In  practice,  we  have 
developed  the  VSTS  method  so  that  these  orbits  may  be  selected  efficiently.)  What  is  of  interest  to  us 
here  is  the  energy  range  for  the  particles. 

Figures  2  and  3  depict  the  energy  diagrams  for  repulsive  and  attractive  potentials.  The  horizontal  lines 
show  typical  energy  bins  for  tracking,  which  are  set  to  span,  with  sufficient  accuracy,  the  Mach  energy 
M^  at  infinity.  The  kinematic  constraint  that  the  kinetic  energy  k  =  e  -  greater  than  zero  sets 
the  integration  limits  on  total  energy: 

repulsive  potential  (x>0)  •  0  <  e  <  oo 
attractive  potential  (x<0)  X  <  E  ^  oo 


Although  the  total  encD  v  range  ;f  to  0  is  kinematically  allowed  for  attractive  potentials,  it  is  normallv 
excluded  from  the  MACH  calculations  because  such  orbits  cannot  connect  to  infinity,  i.e.  they  are 
energetically  trapped.  The  only  way  that  these  orbits  can  become  populated  is  by  collisions.  The 
trapped  electron  model  fills  in  this  kinematic  region  with  an  approximate  equation,  which  is  described 
in  the  next  section. 
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3.  THE  TRAPPED  ELECTRON  MODEL 


It  is  convenient  to  make  a  transformation  from  e  to  kinetic  energy  k  as  the  integration  variable,  then 
the  density  integral  becomes 


where  the  integration  limits  on  k  are 


repulsion'. 

0 

< 

k 

< 

00 

attraction: 

0 

< 

k 

< 

w  for 

trapped  particles 

w 

< 

k 

< 

00  for 

Vlasov  particles 

where  w  =  -x  >  0  is  the  dimensionless  well  depth.  'ITie  variable  k  is  also  indicated  in  Figure  2. 

For  attractive  potentials  (x  <  0)  we  can  build  a  trapped  particle  model  for  the  kinetic  energy 
range  0  <  it  <  ty  by  making  two  assumptions: 

1)  'ITie  source  is  isotropic  (M  =  0).  Hence.  X  =  e  -  k  +  x  Jind  the  density  integral  becomes 

n  =  .iit'^exp(-x)  J|Jexp(- it)it^^(fit(iQ 


2)  All  particles  escape.  Ilius,  we  can  integrate  over  solid  angle  so  that  we  have  J  J  t/Q  =  4k. 
Putting  these  two  assumptions  together  gives 

=  exp(-x)  •  a(«;) 


where  we  define 


a{w)  =  2n~'^^^exp{-k)k^^dk 

In  the  above  equations,  we  have  implicitly  assumed  that  the  temperature  of  the  trapped  and  uiitrapped 
populations  are  the  same.  'Hiis  is  an  approximation,  since  the  higher  energy  particles  are  preferentially 
scattered  out  of  the  well,  leaving  a  lower  temperature  for  the  trapped  component. 

I’o  illustrate  our  normalization,  consider  for  a  moment  the  case  of  repulsive  potentials  (x  >  0).  llie 
only  change  would  be  that  the  limits  on  the  k  integration  would  be  0  to  x  and  we  would  have 
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a(oo)  =  -  T  - 

[2 


=  1  . 


so  that  we  would  recover  the  Boltzmann  equation  for  electron  density. 

We  can  also,  in  this  manner,  make  a  crude  approximation  for  the  Vlasov  electrons  in  the 
range  w  <  k  <  oo-  We  here  keep  assumption  1)  above,  but  weaken  2)  to:  the  solid  angle  for  escape 
is  independent  of  t  This  gives 

rV  =  exp(-x)  •  P(w)  • 

P(ti>)  =  f "exp(-ib)ifc'^dib 

The  quantity  is  the  solid  angle  normalized  to  ./ir 

g^  =  JJdQ/(4it) 

which  is  also  equivalent,  for  an  isotropic  source,  to  the  neutral  particle  density.  This  type  of 
approximation  works  best,  of  course,  for  straight  line  tracks,  i.e.  for  high  velocities  and  weak  electric 
fields.  Since  a(w)  +  P(ti/)  =  1  ,  we  can  write  in  the  form 

=  exp(-x)  •  (1  -  a(u;))  •  g^ 

Iliis  expression  is  very  approximate.  However,  in  our  MIP  simulations,  it  makes  only  a  small 
contribution  over  most  of  the  chamber  (since  g^  drops  off  as  1/r^)  and  we  have  used  it  because  it  avoids 
the  large  CPU  time  required  to  track  the  Vlasov  electrons. 

In  the  above  treatment,  the  electron  density  is  normalized  to  1.  at  the  source,  which  is  the  same  as  for 
the  ion  density.  'Hiis  is  the  correct  boundary  condition  for  space  conditions,  but  is  not  a  good 
approximation  for  our  plasma  chamber  simulations.  There,  we  have  an  ion  thruster  which  produces  a 
plasma  stream  at  the  front  of  the  chamber,  and  a  hot  wire  filament  located  just  in  front  of  the  ion  gun 
which  provides  isotropic  electrons  by  emitting  an  approximately  equal  current.  Thus,  we  have  used  a 
current  balance  condition  to  set  the  normalization  of  the  electron  density.  We  write  the  total  electron 
density  as  n  =  where 
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•  exp(-x)  • 


T\,  =  <^  •  exp(-x)  •  (1  -  a{w))  ■ 

Here,  and  are  separate  normalization  constants.  'Hie  value  of  relative  to  the  ions  is  determined 
by  approximate  flux  balance  at  the  source  region: 

Flux^  =  •  Vj  =  Flux^  =  1/4  •  Vg 

and  therefore 

tv  =  NJN^  ^  4vjv^  = 

wliere  Mj  =  v^Vj  is  the  ion  Mach  number,  T^,  are  temperatures,  m^,  are  masses,  and  v^,  v-  are 
thermal  velocities  (equal  to  (2k'r/m)’'^)  and  =  2n~^^v^  denotes  average  velocity.  Assuming  the 
potential  is  zero  at  the  source,  the  constant  is  undetermined  by  this  condition,  since  a(0)  -  0.  'Hiis 
constant  is  thus  a  free  parameter  in  our  model.  Since  the  trapped  electron  source  is  indirectly  that  of 
the  neutralizing  filament,  we  expect  that  would  be  some  fraction  of  c^. 

4.  NUMERICAL  RESULTS 


In  order  to  show  the  effects  of  the  model,  we  have  integrated  the  function  a(w)  for  various  well  depths. 
Hiis  trapping  function  a(w)  is  given  in  Figure  4.  We  can  see  that  the  function  rises  rapidly  from  zero 
to  one  in  about  6  kT^.  In  Figure  5,  the  trapped  electron  density  calculated  assuming  unit 
normalization  is  shown.  Hiis  curve  shows  an  even  faster  rise  with  respect  to  w  due  to  the  exponential 
factor  exp(w).  It  is  apparent  that  it  only  takes  a  potential  well  of  a  few  kT^  depth  to  produce  electron 
densities  comparable  to  that  of  the  ions.  I'his  is  the  mechanism  by  which  positive  space  charge  at  the 
throat  of  the  chamber  is  effecuvely  neutralized.  Without  this  strong  negative  trapped  orbit  contribution 
to  the  density,  one  tends  to  get  an  unrealistic  and  unstable  build  up  of  positive  charge. 

To  illustrate  the  numerics,  we  show  below  some  results  for  the  MIP  simulation  with  -5KV  probe.  The 
measured  electrostatic  potentials  are  shown  in  Figure  6.  The  spherical  probe  is  situated  at  Z  =  1.0,  R 
=  0.,  behind  an  uncharged  front  disk  centered  aX  Z  =  0.  We  see  that  experimentally,  the  potential  drops 
off  from  high  negative  values  near  the  probe  to  zero,  and  then  goes  positive.  We  interpret  these  positive 
potentials  as  due  to  the  trapped  electrons  in  balance  with  the  streaming  ions.  'Hie  potential  is 
approximately  +4  Volts  at  the  upstream  side  of  the  disk.  'Hie  4  Volt  contour  follows  the  drift  direction 
downstream  (toward  positive  Z)  on  both  sides  of  the  probe.  For  the  MIP  experiment  kT^  =  5  eV,  and 
this  contour  corresponds  to  a  dimensionless  well  depth  of  about  w  -  -x  -  4/5.  From  Figure  4,  we 
see  that  this  gives  a  trapping  fraction  of  -0..\5,  and  from  Figure  5,  an  un-normalized  electron  densi^  of 
-0.75. 
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TRAPPED  ELECTRON  MODEL 


Figure  5.  Trapped  electron  density  alpha{w)  •  exp(w) 


For  Argon  (40  AMU)  flowing  at  approximately  Mach  number  3,  the  current  balance  condition  yields 


=  2ii*'^*3*2''^/(40*1832)’^  -  0.06 

ITie  normalized  trapped  electron  contribution,  assuming  would  be 

n,  -  0.06  •  0.75  -  0.05 


This  is  of  the  right  order  of  magnitude,  but  as  expected,  is  found  to  be  too  high  to  match  the  ion  density 
at  this  p)oint.  To  fine  tune  the  problem,  we  have  taken  the  value  of  the  free  parameter  to  be  less  than 
c^,  the  approximate  percentage  being  determined  by  the  best  fit  to  the  measured  data.  Figures  7,  8,  9, 
and  10  show  the  Poisson  solutions  with  fixed  ion  density  and  a  floating  trapped  electron  term,  and  with 
normalization  set  at  20,  30,  40,  and  50%  of  c^.  These  Figures  indicate  the  sensitivity  of  the  solutions  to 
the  trapped  electron  normalization  constant.  As  the  normalization  factor  is  increased,  a  shallower 
potential  well  is  required  to  bring  in  the  same  negative  charge,  and  the  positive  potentials  thereby  tend 
to  decrease,  as  can  be  seen  in  Figures  7  to  10.  It  appears  that  a  value  of  30%  gives  the  best  match  to 
experimental  values  around  the  region  of  the  probe,  although  a  precise  comparison  cannot  be  made  due 
to  the  observed  asymmetry  in  measured  potentials.  Since  MACH  is  an  axisymmetric  code,  this 
asymmetry  does  not  appear  in  the  simulations. 


5.  CONCLUSIONS 


In  our  MACH  simulations  of  the  MIP  experiment,  we  have  found  that  the  neglect  of  electron  collisional 
effects  leads  to  unphysical  build  up  of  positive  charge  near  the  throat  of  the  chamber.  We  have  described 
a  trapped  electron  model  which  we  have  developed  to  suppress  this  positive  charge.  Physically,  we 
assume  that  chamber  electrons  can  scatter  into  the  potential  well  created  by  the  ion  space  charge.  The 
electrons  thermalize,  with  temperature  equal  to  that  of  the  untrapped  electrons.  The  numerical  effect 
of  the  algorithm  is  as  follows:  Whenever  the  local  potential  goes  positive,  that  region  receives  a  flood 
of  negative  charge  due  to  electrons  that  fall  into  the  potential  well.  This  quenches  any  unphysical  build 
up  of  positive  potential.  The  well  depth  adjusts  so  that  a  self-consistent  Poisson- Vlasov  solution  is 
obtained.  'Hie  model  contains  one  free  parameter  which  is  the  normalization  constant  for  the  trapped 
electron  density.  'Hie  method  requires  the  calculation  of  numerical  quadratures,  but  no  particle  tracking. 
'Hiis  trapped  electron  model  has  beer*  applied  to  the  MACH  simulations  of  the  MIP  chamber 
experiment,  and  enables  us  to  achieve  stable  solutions  with  potential  distributions  that  are  close  to  the 
measured  ones. 
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